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Computing secular motion under slowly rotating quadratic 
perturbation 
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ABSTRACT 

We consider secular perturbations of nearly Keplerian two-body motion under a per- 
turbing potential that can be approximated to sufficient accuracy by expanding it to 
second order in the coordinates. After averaging over time to obtain the secular Hamil- 
tonian, we use angular momentum and eccentricity vectors as elements. The method 
of variation of constants then leads to a set of equations of motion that are simple 
and regular, thus allowing efficient numerical integration. Some possible applications 
are briefly described. 
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1 INTRODUCTION 

If a binary star, or a planet revolving around a star, is 
perturbed by one or more distant relatively slowly moving 
masses the perturbation can be approximated by a tidal 
field. Fast computation of secular perturbations in such a 
potential is thus desirable. Potential applications include the 
Kozai resonance, motion of a well isolated binary in a star 
cluster and cometary motion under the galactic tidal field. 

Several authors have considered the computation of 
com etary motion perturbat io ns du e t o the galactic tidal field 
e.g. | Hei sler and Tremainel il986l) . IWiegert an d Tremaind 
l|l999lhlBrasserl (1200 J) . More recently iFouchardl i2004l) sug- 
gested the use of the time averaged secular perturbing func- 
tion and the Lagrangian equations for com puting the evo- 
lution of orbital elements of a comet. While iFouchard et alJ 
( 2005) compared various ways to compute the evolution of 
orbi tal elements (using Keplerian or Delaunay elements). 

iBreiter and Rataiczakl J2005) used the so called vecto- 
rial elements, i.e. the areal velocity vector and the eccentric- 
ity vector. However, the perturbing potential considered in 
all these studies was a simplified one, in most cases just the 
term causing the force perpendicular to the galactic plane 
was included. 

In this paper we consider a general quadratic potential. 
For more mathematical generality we include the possibility 
of a constant force (although such a term does not appear in 
a tidal field) and also a slow rotation of the quadratic per- 
turbing field is allowed. We form the simple regular equa- 
tions of motion using the vectorial elements. 

Finally it is necessary to stress that the main purpose 
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of this paper is to introduce the simple secular equations of 
motion when a two-body system is perturbed by a (at most) 
quadratic perturbing potential. 



2 PERTURBATIONS OF VECTORIAL 
ELEMENTS 

2.1 Perturbing function 

Consider a system in which the Keplerian motion of a par- 
ticle around a central mass m is perturbed by a perturbing 
potential U(r,t). In addition, a slow rotation (with angular 
velocity N) is assumed. The Hamiltonian may be written 



u 1 2 rn 
H — -v — - — r 
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U(r,t), 



(1) 



where m is the mass of the central body, r is the position 
vector with components (an, x%, X3), v is the velocity vector 
with components (vi,l)i,vz) i.e. r = v (which also can be 
identified as the momentum vector) and a is the semi-major 
axis of the Kepler motion. The function U is the perturbing 
function, assumed at most quadratic in the coordinates. If 
the time dependence of U is due to rotation of the poten- 
tial with a constant angular velocity, then one can consider 
the system in the rotating coordinate system in which the 
perturbing function is a constant. Thus, using the suffix no- 
tation, one may write in place of U 



U 



fiXi + -XidjXj + eijkNiXjVk, 



(2) 



where there is no time dependence, fi are components of 
a (possible) constant acceleration vector / 
second derivative matrix: 
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fi 



dU 

OXi 1^.-0 



Gi 



d 2 u 

dxidxj 



(3) 



The last term in is the scalar product of the angular ve- 
locity vector N of the perturbing potential and the angular 
momentum vector of the moving particle, here written using 
the Levi-Civita symbol e^u- 

Since the matrix Gij is symmetric it can be diagonalized 
by selecting a suitable coordinate system. However, there 
may be situations in which this is not convenient. Therefore 
we consider the general non-diagonal case. In the theory pre- 
sented here that does not add any substantial complications. 

2.2 Averaging the Hamiltonian 

The auxiliary quantities needed in averaging the Hamilto- 
nian are: 

(i) The angular momentum vector per unit of mass (ac- 
tually areal speed) of the particle 

J = r x v. (4) 

(ii) The eccentricity vector (also known as the Runge- 
Lenz- or Laplace vector) 



v x J 



or, for our purposes, the vector 
E — Jvnae. 



(5) 



(6) 



is more convenient. 

(iii) The mean values of the position vector components 



< Xi >= 



E, 



(V 



(iv) The time averaged coordinate products can be writ- 
ten 

< >= i j2 ^ - 4-* J ' + %k EiEi ' (8) 

where is the Kronecker S. This result was easy to derive 
for < x\ > . The generalization was then formed by educated 
guess and use of computer algebra for confirmation. 

The secular perturbing function R —< U > takes the form 

R = U < Xi > +-^-Gij < XiXj > +NiJ t . (9) 
4m 



elements 

A a = ^((E*G**)&i-G«) (10) 
and the vector F 

then the averaged perturbing function can be written in the 
simple form 

R = FiEi + ^JiA, 3 J 3 + -EiBijEj + NiJi, (13) 

where the vectors F, N , and matrices A, B are constants 
in any particular orbit. 



2.3 Equations of motion 

The reason for the choice of the elements J and E was that 
one obtains for their components the simple Poisson bracket 
relations 

{Ji, Jj} = tijkJk, {Ei, Ej} — tijkJk, {Ji, Ej} = e ijk E k . (14) 

With these formulae one gets the equations of motion (Bre- 
iter and Ratajczak(2005) and references therein) 

v r _ _ dR „ dR ,„ „> 

J =-{J,i?}= Jx—+Ex— (15) 

E =-{E,R}= Ex^+Jx^, (16) 

and by II. '-it the partial derivatives of R with respect to the 
vectorial elements are 

dR 
dJ 
dR 
dE 



N + A J 
F+ BE, 



(17) 
(18) 



where we have used the vector-matrix notation. In deriv- 
ing the above equations one should note that J and E 
are constants of motion under the two-body Hamiltonian 
H — —m/(2a). Thus a can be considered a (numerical) con- 
stant since {J, a} = {E,a} — 0. The equations of motion 
have the integrals 



J-E = 0, 



J 2 + E 2 



ma, R = constant. 



(19) 



Here we must remark that the form of the perturbing 
function is not unique. This is because of the larger-than- 
necessary number of variables and because one may e.g add 
terms proportional to J 2 + E 2 and/or to J ■ E which would 
not affect the final expressions for the derivatives of the el- 
ements in 11511 and llbl . 

Note that when the rotation term JV • J is included, 
the result for the evolution of the vectorial elements will be 
obtained in the rotating coordinate system. 



3 APPLICATIONS 

Possible applications of the vectorial element equations in- 
clude: 

(i) Constant force, or rotating constant force. This may 
be restricted to toy applications, but was included in our 
treatment for the sake of generality and because it is math- 
ematically simple. 

(ii) Isolated binary in a star cluster. Assuming the per- 
turbing stars are distant and move relatively slowly, one may 
approximate the perturbing potential by 



U 



E 



2\R k \ 



(Rk ■ r) 

\R k \ 2 



(20) 



where R k are the positions of the perturbers relative to the 
centre-of-mass of the binary. Writing R k = (X k i, X k 2, X k 3) 
we get for the components of the matrix G 



G 



(Sij — 3X k iX k j) , 



(21) 



where the matrix elements may be time dependent, but this 
does not invalidate the equations presented above. 
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(iii) Kozai resonance. This term is often used in connec- 
tion of triple systems in which averaging over the outer orbit 
gives the perturbing potential of the form 

[/ = g|(r 2 -3(s-r) 2 ), (22) 

where 7713 is the mass of the third body, 63 is the semi-minor 
axis of the outer ellipse and z is the unit normal vector of 
the orbital plane of the distant body. In this case we have 

m 3 , 



Gij — 



(23) 



where Zk are the components of the vector z. 

(iv) Cometary motion under the galactic tide. In this case 
the tidal potential is usually written 



U 



^y^G kk xl + N-J, 



(24) 



which is possible by choosing the coordinate system such 
that xz axis is perpendicular to the galactic plane and the 
a;i-axis points to the galactic centre. 

(v) The perturbing potential due to a small disk of mass 
can also be expressed easily in terms of the vectorial el- 
ements, although not in quadratic form. The secular per- 
turbing function for this case is often written in terms of the 
Delaunay elements in the form R = e£~ 3 G~ 3 (l - 3G~ 2 H 2 ), 
where e is a (small) constant and L, G, H are the Delau- 
nay elements. Here L — ^Jma is a constant, G — \ J\ and 
H — z ■ J , where z is the unit direction vector of the sym- 
metry axis of the potential. Thus 



R = eL~ 



! (l-3 



(z-J) 



■)• 



(25) 



Consequently it is easy to include the effect of any such 
potential into the equations of motion for the vectorial ele- 
ments. Note that this kind of a term could be used to ap- 
proximate the effect of the planets (averaged over orbits) 
to comets that come close, but not too close, to the Solar 
System. 



4 NUMERICAL ASPECTS 

We ran some numerical test computations and compared 
results against direct numerical integration of the equations 
of motion for the coordinates. The tests used the tidal field of 
the Galaxy in the form U = | Gkkx\+N-J (parameters 
from lFouchard et alJ ll2005h l and a value for the semi-major- 
axis of an Oort Cloud comet of a ^ 10, 000AU. 

The experiments suggest: 

If one needs results for a long time and infrequent out- 
put is allowed, then one of the best numerical integrators 
available is the Burlirsh-Stoer method (actually the code 
known as difsyi, originally written by Burlirsh). However, 
this is efficient only when run with optimally long steps, 
which may be too long for some applications. 

If one needs frequent output, then the method of choice 
is the implicit midpoint method. It is very simple to imple- 
ment and gives quite high precision , especially it conserves 
the quadratic integrals of motion ([Huang and LeimkuhleJ 
Jl997l) . ISanz-Serna and Calvd J1994I) ). 

In this case one finds that the stepsize equal to one 



period (i.e. 10 6 years) gives a clearly satisfactory accuracy: 
the plots of the vectorial elements from both computations 
(secular and coordinate integration) agree well, except for a 
small phase error. 

However, for a » 20, 000 AU the results differ and the 
reliability of a secular theory becomes questionable. 



5 CONCLUSIONS 



We ha ve generalized the equations of lBreiter and Rataiczakl 
( 2005) who considered the simple case of perturbing function 
of the form (733 z 2 . 

The equations for vectorial elements, even in the case 
of a general quadratic potential are simple and regular, con- 
trary to w hat one obtains when u sing Keplerian or Delaunay 
elements iFouchard et alJ feOOSTO . 

The integration of the secular theory equations is typi- 
cally faster than direct coordinate integration by two orders 
of magnitude. 

The implicit midpoint method may be even faster than 
Burlirsh-Stoer, but there is no simple way of knowing accu- 
racy and optimal stepsize. 

The most reliable integration method is Burlirsh-Stoer 
extrapolation. 

For comets the secular theory results are not accurate 
for a >> 20, 000AU, as comparison with direct integration 
shows (simply: perturbation can be too large). 

For an integration over only one period at a time the 
midpoint method and Burlirsh-Stoer are nearly equally fast. 
Thus in this case the midpoint method may be preferable 
due to its simplicity. 
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